################################################
#	Spry and Nunnally (2025)
# Replication Data for "Modern American Dilemma"
# Exploratory Factor Analysis
################################################

# Front Matter ------------------------------------------------------------

rm(list=ls())		# clear objects in memory

setwd("/Users/amberspry/Dropbox/Research/Data/2016 CMPS/Analysis")
# Data imported from 2016 CMPS

source("2026_sprynunnally_source.R")
library(MVN)
library(kableExtra)
library(psych)
library(tidyr)
library(Hmisc)
library(dplyr)
library(lavaan)


# Create Weighting Vectors ------------------------------------------------

# Clean Data
# Remove rows with missing values to run FA
# Create weighting vector

vars.bw           <- na.omit(vars.bw1)
vars.bw.weight    <- sub.bw$weight[complete.cases(vars.bw1)]

vars.black        <- na.omit(vars.black1)
vars.black.weight <- sub.bw$weight[complete.cases(vars.black1)]

vars.white        <- na.omit(vars.white1)
vars.white.weight <- sub.bw$weight[complete.cases(vars.white1)]


# Define Exploratory Factor Analysis Model --------------------------------

# Un-weighted PCA

pca1 <- prcomp(na.omit(vars.bw), scale. = TRUE)                   # Run a PCA to determine number of underlying components
summary(pca1)                                                     # focus on values in "proportion of variance" - looks like 3 components 
screeplot(pca1, npcs = 9, type = "barplot")                       # Scree plot of PCA - again looks like 3 components
title(main = "Principal Component Analysis")

# Weighted PCA

pca1.weighted <- principal(vars.bw,
                           nfactors = 9,
                           rotate = "none",
                           weights = vars.bw.weights.complete)

plot(pca1.weighted$values, type = "b",
     xlab = "Component Number",
     ylab = "Eigenvalue",
     main = "Principal Component Analysis",
     xaxt = "n")                                                  # Suppress default x-axis
axis(1, at = 1:9, labels = 1:9)                                   # Add custom x-axis with all 9 components
abline(h = 1, lty = 2)                                            # Kaiser criterion line

# Based on summary it appears 3 components exist


# Exploratory Factor Analysis with 3 Factors for ALL Respondents ------------------------------

# Weighted Exploratory Factor Analysis using Psych package

fa2.bw.weighted <- fa(vars.bw,
                       nfactors = 3,
                       rotate = "varimax",
                       weights = vars.bw.weight,
                       fm = "ml")                   # maximum likelihood estimation

# View results
print(fa2.bw.weighted, cut = 0.1)     

# Factor loadings
fa2.bw.weighted$loadings

# Proportion of variance explained
fa2.bw.weighted$Vaccounted


# 3 Component Factor Analysis for BLACK respondents -------------------------

# Weighted Exploratory Factor Analysis using Psych (to explore difference in weight treatment)

fa2.black.weighted <- fa(vars.black,
                      nfactors = 3,
                      rotate = "varimax",
                      weights = vars.black.weight,
                      fm = "ml")                   # maximum likelihood estimation

# View results
print(fa2.black.weighted, cut = 0.1)

# Factor loadings
fa2.black.weighted$loadings

# Proportion of variance explained
fa2.black.weighted$Vaccounted


# 3 Component Factor Analysis for WHITE Respondents -----------------------

# Weighted Exploratory Factor Analysis using Psych (to explore difference in weight treatment)

fa2.white.weighted <- fa(vars.white,
                      nfactors = 3,
                      rotate = "varimax",
                      weights = vars.white.weight,
                      fm = "ml")                   # maximum likelihood estimation

# View results
print(fa2.white.weighted, cut = 0.1)

# Factor loadings
fa2.white.weighted$loadings

# Proportion of variance explained
fa2.white.weighted$Vaccounted


# Table of Uniqueness Scores ----------------------------------------------

fa2.bw.weighted$uniquenesses
fa2.black.weighted$uniquenesses
fa2.white.weighted$uniquenesses
